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1.  INTRODUCTION 

In  this  document,  several  different  Fast  Fourier  Transform  (FFT)  algorithms  are 
presented.  Each  FFT  is  tested  with  respect  to  timing,  accuracy,  storage  re¬ 
quirements,  and  other  restrictions.  The  test  results  for  all  the  FFT  codes  are 
then  summarized  and  compared.  Complete  FORTRAN  codes  and  instructions  for 
their  use  accompany  the  discussions. 


2.  THE  FAST  FOURIER  TRANSFORM 

Several  different  definitions  of  the  discrete  Fourier  transform  appear  in  the 

literature.  The  definition  which  we  shall  use  is  as  follows:  let  xq»xj . 

xn_^  be  a  sequence  of  data,  points.  The  discrete  Fourier  transform  (DFT)  of 
Xq,  x^,  ...,  xn_i  is  the  sequence  X^,X^,  ...,  X^_^  given  by 


2ni1k 


X^-S  xe  n  (k=0,l,  ...,n-l), 
j=0 


where  i=v^l.  The  inverse  discrete  Fourier  transform  (IDFT)  of  X^.X^,  ...,  Xn_^ 


n-1  2iri1k 


j  ^  Xfce  n  (j=0,l . n-1) 


The  algorithm  for  efficiently  computing  the  DFT  and  the  IDFT  is  called  the 
Fast  Fourier  Transform  and  was  rediscovered  in  1965  by  Cooley  and  Tukey  [1]. 
Several  different  FORTRAN  codes  for  computing  the  FFT  will  be  presented.  In 
some  cases,  the  codes  allowed  the  calculation  only  of  the  DFT  and  not  of  the 
IDFT  (or,  in  those  instances  where  the  definition  of  the  DFT  differed  from  ours, 
the  calculation  only  of  the  IDFT) .  The  codes  have  been  modified  so  that  both 
the  DFT  and  the  IDFT  can  be  calculated  by  the  same  subroutine. 
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The  following  data  are  presented  for  the  various  FFT  algorithms: 

(i)  Reference  to  the  literature. 

(ii)  FORTRAN  program  listing. 

(iii)  Timing:  the  average  computation  time  for  an  FFT  of  n  data  points 
on  the  Raytheon  704  computer  where  it  is  assumed  that  n  is  a  power 
of  2. 

(iv)  Accuracy:  if  Xq,x^,  . ..,  xn_^  is  the  input  sequence,  we  compute 
■  Vj  =  IDFt|dFT[x^.]|  -  x^ 

for  j=0,l,  ...,  n  and  hen  calculate  the  root -mean-square  (RMS) 
error 


e  will  be  computed  for  the  sequences 
x^  =  exp  (27iijk/n) 

for  k=2  and  n-64,  128,  256,  512  &  1024. 

(v)  Number  of  executable  FORTRAN  statements 

(vi)  Internal  array  storage  requirements. 

(vii)  E:  rerna1  array  storage  requirements. 

(viii)  Tjpe  of  arithmetic  used:  real,  complex,  or  integer, 
(ix)  Other  restrictions 
(x)  Program  calling  sequence 


2.1  CHARACTERISTICS  OF  THE  FFT  ALGORITHMS 

FFT  Algorithm  I 
(i)  Reference:  [2] 

(ii)  FORTRAN  program  listing:  See  Figure  1. 


( 
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SUBROUTINE  FFTC A#M# N* IS) 

DIMENSION  A(N> 

COMPLEX  A#U*W#T 

N-2**M 

NV2-N/2 

NM1-N-1 

J-l 

DO  30  I  =  1  .*  NM 1 
IF  Cl.  GE.  J>  GO  TO  10 
T-AC  J> 

AC  J)«A<  I  > 

AC  I  >-T 
10  K-NV2 

20  IF  CK.GE.J)  GO  TO  30 
J-J-K 
K-K/2 
GO  TO  20 
30  J-J+K 

PI  -3.14159265 
DO  80  L-1#M 
LE«2**L 
LEl-LE/2 
U- CMPLX  C  1 .  *  0  *  ) 

IF  CIS)  40#  50#  50 

40  V- CMPLX C  COSC  PI/LEl)#-SINCPI /LEI ) ) 

GO  TO  60 

50  W- CMPLX C COSCPI /LEI >#  SINCPI /LEI ) ) 
60  DO  80  J»l#LEl 
DO  70  I-J#N#LE 
IP-I+LE1 
T-AC  IP)*U 
AC IP)-AC I >-T 
70  ACD-ACD+T 
80  U-U*U 
RETURN 
END 

Figure  1.  FFT  Algorithm  I  Program  Listing 


5  January  1972 


4 


System  Development  Corporation 
TM-485  7 / 100/00 


(ill) 

(iv) 


(v) 

(vi) 

(vii) 

(viii) 

(ix) 

(x) 


(i) 

(ii) 

(iii) 

(iv) 


Timing:  .0033n  log2n 
Accuracy: 

Number  of  Data  Points 
64 
128 
256 
512 
1024 


RMS  error 
1.9675  E-06 
3.2880  E-06 
6.4563  E-06 
1.1038  E-05 
1.8687  E-05 


Number  of  executable  FORTRAN  statements:  32 
Internal  array  storage  requirements:  None. 

External  array  storage  requirements:  A  complex  array  of  dimension 
n,  where  n  is  the  number  of  data  points. 

Type  of  arithemtic  used:  Complex 

Other  restrictions:  Number  of  data  points  must  be  a  power  of  2. 
Calling  sequence:  CALL  FFT  (A.MjNjIS),  where 

A  «*  input  array  of  samples 
N  -  2**M  ■  number  of  samples 
M  -  log2(N) 

IS  ■  -1,  forward  transform 

■  +1,  inverse  transform  (unnormalized). 


FFT  Algorithm  II 


Reference:  [3] 

FORTRAN  program  listing 
Timing:  .0086n  log2n. 
Accuracy 

Number  of  Data  Points 
64 
128 
256 
512 
1024 


See  Figure  2. 


RMS  Error 
9.2160  E-07 
1.1040  E-06 
1.2825  E-06 
1.4254  E-06 
1.5463  E-06 


on ooooooo 
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SUBROUTINE  FFT(X«NSTAGE«SI GN) 

INPUT* 

X<  2*1024)  t  DATA  INPUT  IN  COLUMN  I 
NSTAGEt  POWER  OF  TWO  WHICH  N  IS* 

N  ■  2**NSTAGE 

SIGN*  «-l>  FORWARD  TRANSFORM 
«+l>  INVERSE  TRANSFORM 

OUTPUT* 

X<2*1024)*  FORWARD* TRANSFORMED  DATA  OUTPUT  IN  COLUMN  1 
INVERSE* TRANSFORMED  DATA  OUTPUT  IN  COLUMN  2 
COMPLEX  XC  2* 1024)*  W 
INTEGER  R*  SIGN 
N“2**NSTAGE 
N2“N/2 
FLTN-N 

PHI2N-6. 28318 53/FLTN 
DO  3  J- 1>NSTAGE 
N2J-N/C  2**J> 

NR=N2J 
NI«C2**J)/2 
DO  2  I*»1*NI 
IN2J»<  7.-1  )*N2J 
FLIN2J-IN2J 
XSIGN-SIGN 

TEMP  »  FLIN2J*PHI 2M*XSI GN 
W-CMPLXC  COS<  TEMP) *  SI N( TEMP  >  > 

DO  2  R»1*NR 
ISUB-R+IN2J 
ISUB1«R+IN2J*2 
ISUB2-ISUB1+N2J 
I SUB3* I SUB+N2 

X( 2* ISUB)"X( 1* ISUBl >  +  W*X( 1> I SUB2) 

X ( 2> I  SUB 3 ) >X ( 1 > I  SUB 1 )  *  W*X< 1> ISUB2) 

2  CONTINUE 
DO  3  R>LN 

3  XC1*R)®XC2>R> 

IF  (SIGN*LT*0« )  RETURN 
DO  4  h»l>N 

4  X(2*R)*X(  1*R)/FLTN 
RETURN 

END 


Figure  2.  FFT  Algorithm  II  Program  Listing 
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(v) 

(vi) 

(vii) 

(viii) 

(ix) 


(i) 

(ii) 

(iii) 

(iv) 


(v) 

(vi) 

(vii) 


(viii) 


Number  of  executable  FORTRAN  statements:  28 
Internal  array  storage  requirements:  None 

External  array  storage  requirements:  A  2xn  complex  array,  where  n 
is  the  number  of  data  points. 

Type  of  arithmetic  used:  Complex. 

Calling  sequence:  CALL  FFT  C.i,  NSTAGE,  SIGN),  where 
X(2,1024)  ■  data  input  in  column  1 
NSTAGE  =  log2(N),  where  N  =  the  number  of  data  points 
SIGN  =  -1,  forward  transform:  data  output  in  column  1 

=  +1,  inverse  transform:  normalized  data  output  in 
column  2 

FFT  Algorithm  III 
Reference:  [4] 

FORTRAN  program  listing:  See  Figure  3. 

Timing :  . 0026n  log^n . 

Accuracy: 

Number  of  Data  Points  RMS  er^r 

64  1.7072  E-06 

128  1.4605  E-06 

256  2.0886  E-06 

512  2.0916  E-06 

1024  2.3507  E-06 

Number  of  executable  FORTRAN  statements:  478 
Internal  array  storage  requirements:  312 

External  array  storage  requirements:  Two  real  arrays  of  dimension 
n,  representing  the  real  and  imaginary  parts  of  the  input  sequence 
of  length  n. 

Type  of  arithmetic  used:  Real 
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SUBROUTINE  FFTC A*  B*  NfOT*N*NSPAN* I SN> 

C  MULTIVARIATE  COMPLEX  FOURIER  TRANSFORM*  COMPUTED  IN  PLACE 
C  USING  MIXED- RADIX  FAST  FOURIER  TRANSFORM  ALGORITHM. 

C  BY  R.  C.  SINGLETON*  STANFORD  RESEARCH  INSTITUTE*  OCT.  1968 
C  ARRAYS  A  AND  B  ORIGINALLY  HOLD  THE  REAL  AND  IMAGINARY 
C  COMPONENTS  OF  THE  DATA*  AND  RETURN  THE  REAL  AND 
C  IMAGINARY  COMPONENTS  OF  THE  RESULTING  FOURIER  COFFICIENTS. 

C  MULTIVARIATE  DATA  IS  INDEXED  ACCORDING  TO  THE  FORTRAN 
C  ARRAY  ELEMENT  SUCCESSOR  FUNCTION*  WITHOUT  LIMIT 

C  ON  THE  NUMBER  CF  IMPLIED  MULTIPLE  SUBSCRIPTS. 

C  THE  SUBROUTINE  IS  CALLED  ONCE  FOR  EACH  VARIATE* 

C  THE  CALLS  FOR  A  MULTIVARIATE  TRANSFORM  MAY  BE  IN  ANY  ORDER. 
C  NTOT  IS  THE  TCTAL  NUMBER  OF  COMPLEX  DATA  VALUES. 

C  N  IS  THE  DIMENSION  OF  THE  CURRENT  VARIABLE* 

C  NSPAN/N  IS  THE  SPACING  OF  CONSECUTIVE  DATA  VALUES 
C  WHILE  INDEXING  THE  CURRENT  VARIABLE. 

C  THE  SIGN  OF  ISN  DETERMINES  THE  SIGN  OF  THE  COMPLEX 
C  EXPONENTIAL*  AND  THE  MAGNITUDE  OF  ISN  IS  NORMALLY  ONE* 

C  A  TRI-VARIATE  TRANSFORM  WITH  ACN1*N2*N3)*  B(N1*N2*N3) 

C  IS  COMPUTED  BY 

C  CALL  FFTC  A*B*N1*N2*N3*N1*N1*  1) 

C  CALL  FFT<A*B*N1*N2*N3*N2*N1*N2*1> 

C  CALL  FFTC  A*B*N1*N2*N3*N3*N1*N2*N3*  1  > 

C  FOR  A  SINGLE- VARIATE  TRANSFORM* 

C  NTOT  **  N  ■  NSPAN  »  (NUMBER  OF  COMPLEX  DATA  VALUES)*  F.G. 

C  CALL  FFTCA*B*N*N*N*  1) 

C  THE  DATA  MAY  ALTERNATIVELY  BE  STORED  IN  A  SINGLE  COMPLEX 
C  ARRAY  A*  THEN  THE  MAGNITUDE  OF  ISN  CHANGED  TO  TWO  TO 
C  GIVE  THE  CORRECT  INDEXING  INCREMENT  AND  AC 2)  USED  TO 

C  PASS  THE  INITIAL  ADDRESS  FOR  THE  SEQUENCE  OF  IMAGINARY 

C  VALUES*  E.  G. 

C  CALL  FFTC A* AC 2) *NTOT*N*NSPAN*  2) 

C  ARRAYS  ATCMAXF)*  CKCMAXF)*  BTCMAXF)*  SKCMAXF)*  AND  MPCMAXP) 

C  ARE  USED  FOR  TEMPORARY  STORAGE*  IF  THE  AVAILABLE  STORAGE 
C  IS  INSUFFICIENT*  THE  PROGRAM  IS  TERMINATED  BY  A  STOP. 

C  MAXF  MUST  BE  . GE«  THE  MAXIMUM  PRIME  FACTOR  OF  N. 

C  MAXP  MUST  BE  «GT.  THE  NUMBER  OF  PRIME  FACTORS  OF  N. 

C  IN  ADDITION*  IF  THE  SQUARE-FREE  PORTION  K  OF  N  HAS  TWO  OR 

C  MORE  PRIME  FACTORS*  THEN  MAXP  MUST  BE  • GE.  K-l* 

DIMENSION  ACN>*  BCN) 

C  ARRAY  STORAGE  IN  NFAC  FOR  A  MAXIMUM  OF  11  FACTORS  OF  N. 

C  IF  N  HAS  MORE  THAN  ONE  SQUARE- FREE  FACTOR*  THE  PRODUCT  OF  THE 
C  SQUARE- FREE  FACTORS  MUST  BE  *LE*  210 
DIMENSION  NFACC  ID* NPC 209 ) 


Figure  3.  FFT  Algoiltnm  III 
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C  ARRAY  STORAGE  FOR  MAXIMUM  PRIME  FACTOR  OF  83 
DIMENSION  AT< 23)# CK<  23>*  BTC 23>*  SKC  23) 

EQUIVALENCE  CI*II> 

C  THE  FOLLOWING  TWO  CONSTANTS  SHOULD  AGREE  WITH  THE  ARRAf  DIMENSIONS* 
MAXF-23 
MAXP-209 

IFCN  .LT.  2)  RETURN 
INC-ISN 

RAD- 6. 0 * AT AN Cl. 0) 

S72- RAD/5.0 
C72  ■  COS  CS72) 

S72-SINCS78) 

S120-SQRTC0. 75) 

IFCISN  .GE.  0)  GO  TO  10 
S72— S72 
SI  80— SI  20 
RAD— RAD 
INC— INC 
10  NT-INC*NTOT 
KS-INC*NSPAN 
KSPAN-KS 
NN-NT-INC 
JC-KS/N 

RAD F- RAD*  FLO ATC JC)*0*5 
1-0 
JF-0 

C  DETERMINE  THE  FACTORS  OF  N 
M-0 
K-N 

GO  TO  20 
15  M-M+l 

NFACCM)-4 

K-K/16 

20  IF  CK-CK/16)*16«EQ*0>  GO  TO  15 
J-3 
JJ-9 

GO  TO  30 
25  M-M+l 

NFACCM)- J 
K-K/JJ 


(Cont’d)  Figure  3.  FFT  Algorithm  III 
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30  irCMOMlOJJ)  •IQ*  0)  GO  TO  85 
J»J*2 
JJ- J**2 

I  F<  JJ  .LE.  K)  GO  TO  30 
IF<K  • GT.  4)  GO  TO  40 
KT-M 

NFAC<M+  1  )-K 
IF<K  .NE.  1)  M-M+l 
GO  TO  80 

40  IF<K-<K/4>*4  .NE.  0)  GO  TO  50 
M-M+l 
NFAC<M>«2 
K-K/4 
50  KT-M 
J-2 

60  IFCMOIXKjJ)  .NE.  0)  GO  TO  70 
M-M+l 
NFACCM>-J 
K-K/J 

70  J«<  < J+l )/2>*2+l 

I  F<  J  «LE.  K>  GO  TO  60 
80  IF<KT  .EQ.  0)  GO  TO  100 
J-KT 
90  M-M+l 

NFACCM)-NFACCJ> 

J-J-l 

I  F<  J  .NE.  0)  GO  TO  90 
C  COMPUTE  FOURIER  TRANSFORM 
100  SD-RADF/FLOATC  KSPAN ) 

CD-2.0*SIN<SD>**2 

SD-SIN<SD+SD> 

KK-1 
I-I  +  l 

I FTNFACC  I  >  *NE.  2)  GO  TO  400 

C  TRANSFORM  FOR  FACTOR  OF  2  ( INCLUDING  ROTATION  FACTOR) 
KSPAN-KSPAN/2 
Kl-KSPAN+2 
210  K2-KK+KSPAN 
AK-ACK2) 

BK«B< K2) 

A<K2>-A<KK)-AK 

B<K2)-B<KK)“BK 

A<KK)-A<KK)+AK 

B<KK)-BCKK>+BK 


(Cont'd)  Figure  3.  FFT  Algorithm  III 
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KK-1C2+KSPAN 

IFCKX  *LE.  NN)  GO  TO  210 
KK»KK-NN 

IFCKK  »LE«  JC)  GO  TO  210 
IFCKK  *GT*  KSPAN)  GO  TO  800 
220  C1-1.0-CD 
Sl-SD 

230  K2-KK+KSPAN 

AK”ACKK)-ACK2) 

BK-BCKK)-BCK2) 

ACKK)«ACKK)+ACK2) 

BCKK)«BCKK)+BC K2) 

ACK2)-C1*AK-S1*BK 

BCK2)«S1*AK+C1*BK 

KK«K2*KS?AH 

IFCKK  *LT*  NT)  GO  TO  230 

K2-KK-NT 

Cl  —  C2 

KK-K1-K2 

IFCKK  e GT.  K8)  GO  TO  230 
AK-Cl-CCr*Cl+SD*Sl) 

S1BC  SD*C1-CD*S1 )  +  Sl 

C  THE  FOLLOWING  THREE  STATEMENTS  COMPENSATE  FOR  TRUNCATION 
C  ERROR.  IF  ROUNDED  ARITHMETIC  IS  USED*  SUBSTITUTE 

C  Cl 3 AX 

C1«0. 5/C AK**2+Sl**2)+0. 5 

S1-C1*S1 

C1-C1*AK 

KK-KK+JC 

IFCKK  .LT.  K2)  GO  TO  230 

K1-K1+INC+INC 

KK«  C  K 1 -KSPAN) / 2+ JC 

IF  CKK.LE* JC+JC)  GO  TO  220 

GO  TO  100 

C  TRANSFORM  FOR  FACTOR  OF  3  C OPTIONAL  CODE) 

320  Kl-KK+KSPAN 
K2-K1+KSPAN 
AK-ACKK) 

BK-BCKK) 

AJ"  AC  K1  )+AC  K2) 

BJ«BC  Kl  )+BCK2) 

AC  KK) “AK+AJ 
BCKK)»BK+BJ 
AK— 0* 5*AJ+AK 
BK—  0.5*BJ*BK 


(Cont'd)  Figure  3.  FFT  Algorithm  III 
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AJ-CACKD  -ACK2)  )  +  S120 

BJ»(B(K1 )-fl(K2> )*S120 

A<  K1 ) ■ AK-B J 

BCKD-BK+AJ 

A<K2) “AK+BJ 

B<K2)»BK-AJ 

KX-K2+KSPAN 

IFCKK  .LT.  NN>  GO  TO  320 
KK*»KK-NN 

IFCKK  .LE.  KSPAN)  GO  TO  320 
GO  TO  700 

C  TRANSFORM  FOR  FACTOR  OF  4 

400  IF(NFACCI)  .NE.  4>  GO  TO  600 
KSPNN=KSPAN 
KSPAN-KSPAN/4 
410  Cl*1.0 
Sl-0.0 

420  Kl-KK+KSPAN 
K2-K1+KSPAN 
K3-K2+KSPAN 
AKP-ACKK)+A<K2) 
AXM«A(KK)-ACK2) 
AJP-A<K1)+A<K3) 
AJM-A<K1)-A<K3) 

A( KK ) "AKP+AJP 
AJP-AKP-AJP 
BKP*B(KK)+B(K2) 
BKM«B<KK>-B<K2> 

BJP»B<Ki )+B< K3> 
BJM-B<Ki>-BCK3) 

B(KK>*BKP+BJP 

BJP-BKP-BJP 

1FCZSN  .LT.  0>  GO  TO  450 

AKP-AKM-BJM 

AKM-AKM+BJM 

BKP-BXM+AJM 

BKM-BKM-AJM 

XFCS1  .EQ.  0.0)  GO  TO  460 
430  A<KI)-AKP*C1-BKP*S1 
B<K1)»AKP*S1+BKP*C1 
A<K2)»AJP*C2-BJP*S2 
B(K2> -AJP*S2+BJP* C2 
A<K3)-AKM*C3-BKM*S3 
B(K3)«AKM*S3+BKM*C3 
KK-K3+XSPAN 

IFCKK  .LI.  NT)  CO  TO  420 

(Cont'd)  Figure  3.  FFT  Algorithm  III 
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/ 


440  C2«C1-<C1»C1+SD*S1) 

S1-(SD*C1-CD*S1)+S1 

C  THE  FOLLOWING  THREE  STATEMENTS  COMPENSATE  FOR  TRUNCATION 
C  ERROR*  IF  ROUNDED  ARITHMETIC  IS  USED#  SUBSTITUTE 

C  C1-C2 

Cl-0.5/(C2**2+Sl**2>+0*5 

51- C1*S1 
C1-C1*C£ 

C2«C1**2-S1**2 

S2>2«0*C1*S1 

C3«C2*C1-S2*S1 

S3-C2*S1**S2*C1 

KK-KK-NT+JC 

IFCKK  *LE*  KSPAN  >  GO  TO  420 
KK-KK-KSPAN+INC 
I  F(KK  *LE*  JC)  GO  TO  410 
I FC KSPAN  *EQ*  JC>  GO  TO  800 
GO  TO  100 
450  AKP-AKM+BJM 
AK4-AKM-BJM 
BKP-BKM-AJM 
BKM-BKM+AJM 

I  Ft  SI  *NE*  0*0)  GO  TO  430 
460  ACKD-AKP 
BCK1)»BKP 
ACK2)-AJP 
BCK2)-BJP 
ACK31-AKM 
B(K3)-BKM 
KK-K3+KSPAN 

IFCKK  *LE*  NT)  GO  TO  420 
GO  TO  440 

C  TRANSFORM  FOR  FACTOR  OF  5  (OPTIONAL  CODE) 

510  C2»C72**2-S72**2 

52- 2*0*C72*S72 
520  Kl-KK+KSPAN 

K2-K1+KSPAN 

K3-K2+KSPAN 

K4-K3+KSPAN 

AKP-ACKD+ACK4) 

AKM-ACK1)-A<K4) 

BKP-BCKU+BCK4) 

BKM  -  BCKD-BCK4) 


(Cont’d)  Figure  3.  FFT  Algorithm  III 
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AJP-ACK1»+ACK3> 

AJM-ACX2>-ACK3) 

BJP-B(K2>+B<K3> 

BJM-B(K2>-B<K3> 

AA-A(KK) 

BB-BOOO 

A(  KX  )  «  AA*AKP+AJP 

BCKK)  -BB+BKP' BJP 

AK-AKP*C72+AJP*C2+AA 

BK-BKP*C78+BJP*C2+BB 

AJ-AKM*  S  72+AJM*  58 

BJ-BKM*S72+BJM*S2 

ACK15-AX-BJ 

A(K4)-AK+BJ 

B<K1>»BK+AJ 

B(K4)«BK-AJ 

AK-AKP* C2+AJP* C  78+ AA 

BK-BKP*  C2+BJP*  C78+BB 

AJ-AKM* S2-AJM*S72 

BJ-BKM* S2- BJM*  STB 

A(K2>-AK-BJ 

A(K3)-AK+BJ 

B(K2)-BK>AJ 

BCK3>-BK-AJ 

KK-K4+KSPAN 

I  F(KK  »LT*  NN>  60  TO  520 
KK-KK-NN 

IFCKK  .LE.  KSPAN)  60  TO  520 
60  TO  700 

C  TRANSFORM  FOR  ODD  FACTORS 
600  K-NFACC I > 

KSPMN-KSPAN 

KSPAN-KSPAN/K 

IF(K  •  EQ*  3)  60  TO  320 

IF(K  *EQ*  5>  60  TO  510 

1F(K  *EQ.  JF)  60  TO  640 

JF-K 

Sl-RAD/FLOATCK) 

Cl-COS(Sl) 

Sl-SZN(Sl) 

IFCJF  •  6T*  MAXF)  60  TO  996 
CKC  JF)-1.0 
SK( JF)-0*0 
J-l 


(Con  :'d)  Figure  3.  FFT  Algorithm  III 
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630  CKC  J>»CKCK>*C1+SK<K>*S1 
SK(  J)-CK(K)*S1-SK<K)*C1 
K-K-l 

CKC K> -CKC J> 

SKCK>-«SK( J> 

J-J+l 

IFCJ  .LT.  K*  60  TO  630 
640  Kl-KK 

K2-KK+KSPNN 

AA-A(KK) 

BB-BCKK) 

AK-AA 

BK-BB 

J«1 

K1-K1+KSPAN 
650  KP-K2-KSPAN 
J-J+l 

ATC J>-ACK1>+ACK2> 

AK-ATC J5+AK 
BT< J)-BCK1>+B<K2> 

BK-BTC J)+BK 
J-J+l 

AT(  J  )  -  A(  K 1 )  ■' AC  K2  ) 

BTC  J/-BCK1  )-BCK2> 
K1-K1+KSPAN 

IFCK1  .LT*  K8)  GO  TO  650 

ACKK1-AK 

BOOO-BK 

Kl-KK 

K2-KK+KSPNN 

J»1 

660  K1-K1+KSPAN 
K2-K2-KSPAN 
JJ-J 
AK-AA 
BK-BB 
AJ-0.0 
BJ-0.0 
K-l 


(Cont’d)  Figure  3.  FFT  Algorithm  III 
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670  K-K+l 

AK-AT(K)*CK< JJ)+AK 
BK-BT(K)*CK< JJ)+BK 
K-K+l 

A J- AT< K ) *  SK<  J J ) +AJ 
B  J-  BT(K)*SK<  JJ)+BJ 
JJ* JJ+J 

IF<  JJ  .GT.  JF)  JJ-JJ-JF 

IF(K  .LT.  JF)  GO  TO  670 

K-JF-J 

A(K1)-AK-BJ 

B(K1  )-BK+AJ 

A< K2)«AK+BJ 

B< K2J-BK-AJ 

J-J+l 

I  F< J  .LT*  K)  GO  TO  660 
KK-KK+KSPNN 

I  F(KK  .LE.  NN)  GO  TO  640 
KK-KK-NN 

I  F(KK  »LE*  KSPAN)  GO  TO  640 

C  MULTIPLY  BY  ROTATION  FACTOR  (EXCEPT  FOR  FACTORS  OF  2  AND  4) 
700  I F( I  »EQ.  M)  GO  TO  800 
KK-JC+1 
710  C2-1.0-CD 

51- SD 
720  C1-C2 

52- S1 

KK-XK+KSPAN 
730  AK-A(KK) 

A(  KK) »C2*AK- S2*B< KK) 

B<  KK) -S2+AK+C2* B( KK) 

KK-KK+KSPNN 

I  F(KK  .LE.  NT)  GO  TO  730 

AK-S14S2 

S2-S1*C2+C1*S2 

C2-C14C2-AK 

KK-KK-NT+KSPAN 

IFCKK  .LE.  KSPNN)  GO  TO  730 

C2-C1"(CD*C1+SD*S1) 

Sl-Sl+C  SD*C1-CD+S1 ) 

C  THE  FOLLOWING  THREE  STATEMENTS  COMPENSATE  FOR  TRUNCATION 
C  ERROR.  IF  ROUNDED  ARITHMETIC  IS  USED*  THEY  MAY 
C  BE  DELETED. 


’■> 
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Cl»0.5A'C8**2+Sl**2>+0.5 

S1-C1*SI 

C2-C1*C2 

KK-KK-KSPNN+JC 

IFCKK  .LE.  KSPAN)  GO  TO  720 

KK-KK-KSPAN+JC+INC 

I F<KK  .LE.  JC+JC)  GO  TO  710 

GO  TO  100 

C  PERMUTE  THE  RESULT?  TO  NORMAL  ORDER - DONE  IN  TWO  STAGES 

C  PERMUTATION  FOR  SQUARE  FACTORS  OF  N 
SOO  NP< 1 >«KS 

IFCKT  .EG.  0)  GO  890 
K-KT+KT+1 
I F(M  .LT.  K>  K-K-l 
Jnl 

NPCK+D-JC 

810  NPC J+1)-NP( J)/NFACCU> 

MPC K) -NPC K+ 1 > *NFAC( J) 

J-J+l 

K-K-l 

IFCJ  .LT.  K>  GO  TO  810 
K3-NPCK+1) 

KSPAN- NPC 2) 

KK-JC+1 
K2-KSPAN+ 1 
J-l 

IFCN  .NE.  NTOT)  GO  TO  8  50 

C  PERMUTATION  FOR  SINGLE- VARIATE  TRANSFORM  (OPTIONAL  CODE) 
820  AK-ACKK) 

A(KK>-ACK2) 

AC K2>-AK 
BK-BCKK) 

BCKK)-B(K2) 

BCK2)-BK 

KK-KK+INC 

K2-KSPAN+K2 

IFCK2  .LT.  KS)  GO  TO  820 
830  K2-K2-NPC  J) 

J-J+l 

K2-NPC J+l)+K2 

IFCK2  .GT.  NPC  J> )  GO  TO  830 
J- 1 


(Cont’d)  Figure  3.  FFT  Algorithm  III 
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840  IF(KK  .LT.  K8)  GO  TO  820 
KK-KK+INC 
K2-KSPAN+K2 

IF<K2  .LT.  KS)  GO  TO  840 
I F(KK  .LT.  KS)  GO  TO  830 
JC-K3 
GO  TO  890 

C  PERMUTATION  FOR  MULTIVARIATE  TRANSFORM 
850  K-KK+JC 
860  AK-A(XK) 

A<KK)>A<K2) 

A(K2)-AK 

BK-B(KK) 

BCKK)«B(K2> 

B( K2)-BK 

KK-KK+INC 

K2-K2+INC 

I  F(KK  .LT.  K)  GO  TO  860 

KK-KK+KS-JC 

K2-K2+KS-JC 

I FC XK  .LT*  NT>  GO  TO  850 

K2-K2-NT+KSPAN 

KK-KK-NT+JC 

IFCK2  .LT.  KS)  GO  TO  850 
870  K2«K2-NP<J) 

J-J+l 

K2-NPC  J+ 1 )  +K2 

IF(K2  .GT*  NPC  J) )  GO  TO  8  70 
J-l 

880  IF<KK  *LT.  K2>  GO  TO  8  50 
KK^KK+JC 
K2-KSPAN+K2 

IFCK2  .LT*  KS)  GO  TO  880 
I  F<KK  .LT.  KS)  GO  TO  870 
JC-K3 

890  IF<  2*KT+1  . GE«  M)  RETURN 
KSPNN-NP( KT+ 1 ) 

C  PERMUTATION  FOR  SQUARE- FREE  FACTORS  OF  N 
J-M-KT 
NFACCJ+D-l 


(Cont'd)  Figure  3.  FFT  Algorithm  III 
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900  NFACC  J) -NFACC J)*NFACC J+ 1  > 

J-J-l 

IFCJ  .NE.  KT>  GO  TO  900 

KT-KT+ 1 

NN-NFACC  KT>  - 1 

IFCNN  .GT.  MAXP>  GO  TO  998 

JJ«0 

J-0 

GO  TO  906 

905  JJ-JJ-K2 
K2-KK 
K-K+l 

KK-NFACCK) 

904  JJ-KK+JJ 

IFC  JJ  •  GE.  K2)  GO  TO  902 

NPCJ) «JJ 

906  K2-NFACCKT) 

K-KT+ 1 
KK-NFACCK) 

J-J+l 

IFCJ  .LE.  NN)  GO  TO  904 

C  DETERMINE  THE  PERMUTATION  CYCLES  OF  LENGTH  GREATER  THAN  1 
J-0 

GO  TO  914 
910  K-KK 

KK-NPCK) 

NPCK)  —  KK 

IFCKK  .NE.  J)  GO  TO  910 
K3-KK 
914  J-J+l 

KK-NPCJ) 

IFCKK  .LT.  0)  GO  TO  914 
IFCKK  .NE.  J)  CO  TO  910 
NPCJ)— J 

IFCJ  .NE.  NN>  GO  TO  914 
MAXF-INC*MAXF 

C  REORDER  A  AND  B*  FOLLOWING  THE  PERMUTATION  CYCLES 
GO  TO  950 
924  J-J-l 

IFC  NPCJ)  .LT.  0)  GO  TO  924 
JJ— JC 


iS 


(Cont'd)  Figure  3 
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926  KSPAN-JJ 

IFCJJ  .GT.  MAX  F)  KSP AN-MAX  F 

JJ-JJ-KSPAN 

K-NPC J) 

KK»JC*K+II+JJ 

Kl-KK+KSPAN 

K2-0 

928  K2-K2+ 1 

ATCK2)-ACK1) 

BT<K2)«B<K1> 

K1-K1-INC 

1FCK1  .WE.  KK)  GO  TO  928 
932  Kl-KK+KSPAN 

K2«K1-JC*<K+NP(K)> 

K--NP(K) 

936  ACK1>-A<K2> 

BCKD-BCK2) 

K1-K1-INC 

K2-K2-INC 

IFCK1  .NE.  KK)  GO  TO  936 
KK-K2 

I F<K  .NE.  J)  GO  TO  9  32 

Kl-KK+KSPAN 

K2-0 

940  K2-K2+1 

ACK1 >-ATCK2) 

B<K1)-BT<K2) 

K1-K1-INC 

IFCK1  .NE.  KK)  GO  TO  940 
IF< JJ  .NE.  0)  GO  TO  926 
I F< J  .NE.  1)  GO  TO  924 
950  J-K3+1 

NT-NT-KSPNN 

II-NT-INC+1 

I F<NT  . GE.  0)  GO  TO  924 
RETURN 

C  ERROR  FINISH*  INSUFFICIENT  ARRAY  STORAGE 

998  ISN-O 

C  PRINT  999 
STOP 

999  FORMATC 44H0ARRAY  BOUNDS  EXCEEDED  VI THIN  SUBROUTINE  FFT) 
END 


c 


(Cont  d)  Figure  3.  FFT  Algorithm  III 
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Other  restrictions:  This  algorithm  does  not  require  that  the  number 

n  of  sample  points  be  a  power  of  2;  indeed,  the  prime  factorization 
ki  k 

of  n  =  p^  • • •  Pj  J  need  only  have  the  property  that  p^s  23  for  i=l, 
•  •  •  *  j  • 

Calling  sequence:  CALL  FFT  (A,  B,  N,  N,  N,  ISN),  where 
A  =  array  containing  real  parts  of  input  data, 

B  =  array  cont;  ning  imaginary  parts  of  input  data, 

N  =  number  of  data  points 
ISN  =  -1,  forward  transform 

=  +1,  inverse  transform  (unnormalized). 

For  the  use  of  this  algorithm  for  other  than  radix  2,  see  the  pro¬ 
gram  listing  in  Figure  3. 

FFT  Algorithm  IV 
Reference:  [5] 

FORTRAN  program  listing:  see  Figure  4. 

Timing:  .0042n  log2n 
Accuracy: 

Number  of  Data  Points  RMS  error 

64  9.5800  E-07 

128  1.2362  E-06 

256  1.3819  E-06 

512  1.4222  E-06 

1024  1.4869  E-06 

Number  of  executable  FORTRAN  statements:  18 
Internal  array  storage  requirements:  None 

External  array  storage  requirements:  Two  real  arrays  of  dimension 
n,  representing  the  real  and  imaginary  parts  of  the  input  sequence 
of  length  n. 

Type  of  arithmetic:  Real 
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SUBROUTINE  FFTCX*Y*M,N#IS> 

C  X*Y  -  REAL  &  IMAGINARY  PARTS  OF  THE  INPUT  SEQUENCE 
DIMENSION  XCN>*Y<N> 

DO  10  L0«1*M 
LMX-2**(M-L0> 

LIX«2*LMX 
SCL-6.283185/LIX 
DO  10  LM-1*LMX 
ARG*»<LM“  1  >  *SC". 

C-COS(ARG) 

S— FLOAT<  IS)*SIN<ARG> 

DO  10  LI"*LIX#N«LIX 
Jl-LI-LIX+LM 
J2-J1+LMX 
T1»XCJS)-X<J2) 

T2-YC Jl)-Y< J2) 

X<J1)»XC  J1)*X(  J2> 

Y<J1)*YC  Jll+YC  J2> 

X<J2)-C*T1+S*T2 
10  Y< J2)-C*T2-S*T1 
RETURN 
END 


Figure  4.  FFT  Algorithm  IV  Program  Listing 
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(ix)  Other  restrictions:  Number  of  data  points  must  beea  power  of  2. 

Also,  program  must  be  used  in  conjunction  with  subroutine  RBITS  to 
appropriately  post-order  the  computed  DFT  sequence  (see  Figure  5 
for  a  listing  of  RBITS).* 

(x)  Calling  sequence:  CALL  FFT  (X,  Y,  M,  N,  IS),  where 
X  =  array  containing  real  parts  of  input  data, 

Y  =  array  containing  imaginary  parts  of  input  data, 

N  =  2**M  =  number  of  data  points 
IS  =  -1,  forward  transform 

=  +1,  inverse  transform  (unnormalized) 

The  above  call  to  the  FFT  must  be  followed  by  either 
CALL  RBITS  (X,  Y,  M,  N) 
or  the  two  calls 

CALL  RFLBTS  (X,  N) 

CALL  RFLBTS  (Y,  N) . 


*RBITS  is  designed  to  perform  a  "bit-reflection",  and  not  a  bit-reversal  as 
indicated  by  Markel.  A  simpler  and  more  flexible  bit-reflection  code  is 
shown  in  Figure' 6;  this  program  will  be  recognized  as  the  first  section  of  FFT 
algorithm  i. 
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SUBROUTINE  RBITS<X*Y*M*N) 

C  PERFORMS  XN-PLACE  BIT  REVERSAL  FOR  N«2**M  VALUES  XCI)* 

C  WHERE  M  IS  LESS  THAN  OR  EQUAL  TO  10 

C  OUTPUT  SEQUENCE  IS  Y<I) 

DIMENSION  X<N)*Y<N)*L< 10) 

EQUIVALENCE  <L10*LC  1 )  )*  <L9*LC  2>  )*  <L8*L<3>  )>  <L7*L<  4)  )*  <L6*L<  5)  ) 
EQUIVALENCE  <L5*L<  6))*<L4*L<  7>  )*  CL3*L<8)  )*  CL2*L<9)  )*  CL1*L<  10)) 
DO  20  J-l*10 
L< J)al 

IF  <J-M)  10*  10# 20 
10  L< J)»2**CM+1-J) 

20  CONTINUE 
JN-1 

DO  50  J1»1*L1 
DO  50  J2-J1*L2*L1 
DO  50  J3>J2*L3*L2 
DO  50  J4«J3*L4*L3 
DO  50  J5-J4*L5*L4 
DO  50  J6«J5*L6*L5 
DO  50  J7-J6*L7*L6 
DO  50  J8»J7*L8*L7 
DO  50  J9*J8*L9*L8 
DO  50  JR- J9*L  10*1,9 
IF  CJN-JR)  30*30*40 
30  R-XCJN) 

X<  JN)»XC  JR) 

XC JR)»R 
Fl-YCJN) 

Y< JN)«Yt JR) 

Y< JR)«F1 
40  JN-JN+1 
50  CONTINUE 
RETURN 
END 


Figure  5.  Subroutine  RBITS  Program  Listing 
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SUBROUTINE  RFL3TSCA#N> 

C  PERFORMS  IN-PLACE  'BIT-REFLECTION*  FOR  N-2**M  VALUES  AC  I ) 
C  E*G**  FOR  N-8>  THE  FOLLOWING  MAPPING  WOULD  TAKE  PLACES 
C  0-000  IS  MAPPED  INTO  000 

C  1-001  IS  MAPPED  INTO  100 

C  8-010  IS  MAPPED  INTO  010 

C  3-011  IS  MAPPED  INTO  110 

C  4-100  IS  MAPPED  INTO  001 

C  5-101  IS  MAPPED  INTO  101# ETC* 

C  ORIGINAL  ORDERING  OF  SEQUENCE  IS  DESTROYED 
DIMENSION  ACN) 

NV2-N/8 

NM1-N-1 

J-l 

DO  30  I-1#NM1 

IF  C 1 *GE* J)  GO  TO  10 

T-ACJ) 

AC J)-AC I ) 

ACI)-T 
10  K-NV8 

80  IF  CK.GE.J)  GO  TO  30 
J-J-K 
K-K/2 
GO  TO  80 
30  J-U+K 
RETURN 
END 


Figure  6.  Subroutine  RFLBTS  Program  Listing 
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2.2  FFT  SUMMARY  AND  COMMENTS 

A  summary  of  the  characteristics  of  the  FFT  algorithms  is  given  in  Table  1. 

We  see  that  algorithm  Ill  is  the  fastest  and  most  flexible  with  respect  to  the 
number  of  allowable  data  points  but  requires  more  program  storage  than  the 
others.  The  most  accurate  algorithm  is  II,  but  it  is  about  two  or  three  times 
slower  than  the  others.  Algorithms  I  and  IV  are  the  best  of  the  four  with  re¬ 
spect  to  timing,  accuracy,  and  storage  requirements,  and  since  timing  is 
usually  the  chief  requirement  in  speech  processing,  algorithm  I  appears  to  be 
opt imal . 


Table  1.  .Summary  of  FFT  Characteristics 


5  January  1972 


27 

(last  page) 


System  Development  Corporation 
TM-485 7/100/00 


REFERENCES 


[1]  COOLEY,  J.  W.  and  TUKEY,  J.  W.,  "An  Algorithm  for  the  Machine  Calculation 
of  Complex  Fourier  Series,'1  Math.  Comp. ,  Vol.  19  (1965),  pp.  297-301. 

[2]  COOLEY,  J.  W.  LEWIS,  P.  A.  W.  and  WELCH,  P.  D. ,  "The  Fast  Fourier  Trans¬ 
form  and  its  Applications,"  IEEE  Transactions  on  Education,  Vol.  12  (1969), 
pp.  27-34. 

[3]  UHRICH,  M.  L , ,  "Fast  Fourier  Transforms  Without  Sorting,"  IEEE  Transactions 
on  Audio  and  Electroacoustics,  Vol.  AU-17  (1969),  pp.  170-172 

[4]  SINGLETON,  R.  C.,  "An  Algorithm  for  Computing  the  Mixed  Radix  Fast  Fourier 
Transform,"  IEEE  Transactions  on  Audio  and  Electroacoustics,  Vol.  AU-17 
(1969),  pp.  93-103. 

[5]  MARKEL,  J.  D.,  "FFT  Pruning,"  IEEE  Transactions  on  Audio  and  Electro¬ 
acoustics,  Vol.  AU-19  (1971),  pp.  305-311. 


